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Abstract 

Violation of (semi)-detailed balance conditions in lattice gas automata 
gives rise to unstable spatial fluctuations that lead to phase separation 
and pattern formation in spinodal decomposition, unstable propagating 
modes, driven diffusive systems and unstable uniform flows. 



1 Introduction 



Why is violation of detailed or semi-detailed balance [ Btueckclberg 1952 | of inter- 
est, and what are its consequences? In the context of lattice gas automata viola- 
tion of the Stueckelberg con dition has be en used to increase the Reynolds num- 
ber in hydrodynamic flows [ Hcnon 1992 1 and to model the dynamics of phase 
separation and spatial instabilities |Rothman 1988, Appert 1990, Gerits 1993, 
Alexander 1992 . Lack of detailed balance occurs in several areas of statisti- 



cal and chemical physics, such as driven diffusive systems [2ia 1993|, in Smolu- 
chowski's theory of rapid coagulation [Drake 1972, Ernst 1986| and in the theory 
of neural networks and pattern recognition [ [Clark 1988| , [Perrida 19"90C , where 
master equations with asymmetric transitions rates are being used. 

As soon as detailed balance is lacking, there does not exist an H-theorem. 
In the thermodynamic limit the system may or may not approach a unique 
stable spatially uniform equilibrium state. In the latter case the system seems to 



Bussemaker 1993a |; in the 


; equal times [ 


Henon 1992 



Bussemaker 1992 . Very little is known about the possible stationary states and 
further properties of non-detailed balance models. 

In this paper we study lattice gases that violate the Stueckelberg condition, 
using mean field theory. We are particularly interested in spatial instabilities, 
that lead to phase separation and formation of spatial patterns. The main goal 
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is to show how a Hnear stabiUty analysis or study of the long wavelength hydro- 
dynamic and diffusive modes determines the stability thresholds (temperature, 
range of forces), the initial patterns, i.e. wavelength and direction of maximum 
growth, and the onset time of spatial instabilities. The coarsening of the initial 
patterns, the growth of correlation length and the scaling laws describing the 
late stages of phase separation can not be explained on the basis of this linear 
theory. 

The paper is organized as follows. In section 2 we review the lattice Boltz- 
mann equation, the decay or growth of respectively stable or unstable spatial 
fluctuations, and the Cahn-Hilliard theory. Section 3 presents various mecha- 
nisms for phase separation in lattice gas automata: thermodynamic instabili- 
ties, unstable propagating modes and fluctuating bias fields. In section 4 the 
discussion is extended to instabilities in systems where the spatial symmetry 
is externally broken by an applied field, i.e. driven diffusive systems, or by an 
imposed uniform flow velocity. 



2 Linear stability analysis 
2.1 Lattice Boltzmann equation 

Starting from the nonlinear Boltzmann equation for lattice gas automata with 
local and nonlocal collision rules, we investigate the linear stability of long wave- 
length fluctuations around the spatially uniform state. If the stability criterion 
is not satisfied, long wavelength fluctuations (e.g. concentration fluctuations in 
spinodal decomposition) are unstable and the Cahn-Hilliard theory is used to 
discuss the initial behavior of the structure function and the density-density 
correlation function. In standard notation the lattice Boltzmann equation is 
given by, 

/,(f +cl,t+ 1) - Mr,t) = /,(/(t)) (1) 

where i = 1, 2, . . . , 6 denotes a moving particle state with \ci\ — vq and i = 
denotes a rest particle state with velocity cq = 0. The collision operator Jj 
contains in general not only the standard local collisions, but also nonlocal 
ones. 

We are interested in the stability of long wavelength spatial fluctuations of 



the form |Bussemaker 1993a] , 



SMr, t) = Mr, t) - /f = exp[ifc • f + zx{k)t\ (2) 

around the spatially uniform state /°, given by the solution of Ii{f°) = 0. 
Inserting (|^ into (|^) yields after linearization an eigenvalue equation for the 
A-th mode ijjxj with eigenvalue z\, i.e. 

b 

exp[zx]ipxi{k) = ^exp[~ik ■ Ci]{l + n{k))ijipxj{k) (3) 

3=0 
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where r2(fc) ~ il + ikh. . . .{k 0) is the Fourier transform of the hnearized 
nonlocal coUision operator Ii{f° + Sf). 

What type of information can be obtained from Eq.(H)? We mainly focus 
on the eigenvalue spectrum zx{k) and in particular on the soft hydrodynamic 
modes where z\{k) — > as fc — > 0, that are related to the conservation laws. If 
RezA(fc) > the mode is unstable and the corresponding eigenfunction ipx is 
the order parameter. 

The eigenvalues zx{k) can be determined either numerically or analytically 
by perturbation theory for small |fc| = k. In the long wavelength limit the 
hydrodynamic eigenvalues have the form, 

zx{k) = ikvx + {ikfVx + ... (fc ^ 0) (4) 

In a single component d-dimensional lattice gas with mass and momentum con- 
servation, there are {d + 1) slow modes: two propagating sound modes (labeled 
A = (7 = ±), where = F is the sound damping constant, and Wo- = is 
the speed of sound, defined through — dp/ dp. The remaining [d — 1) shear 
modes (labeled A =_L) are degenerate where Fx = is the kinematic viscosity 
and v_\_ =Q. In a binary mixture To — D \s the diffusion coefficient. 



2.2 Hydrodynamic modes 

We first consider the analytic method. By applying the perturbation theory of 
Gerits et al. [1993] to nonlocal collision operators, ri(fc) = SI + ikK + . . ., the 
damping coefficient of the A-th mode is found to have the form, 

^^^{—- \) ~ -0a|A|5a) (5) 

\LOx 2/ bJx 

or a linear combination of such terms. We have used the notation {a\b) = 
Qihi and {a\M\h) = J^ij o-iMijbj. The current jx is a right eigenvector of the 
collision matrix, fljx = — wja with < ux < 2, and jx is the left eigenvector 
corresponding to the same eigenvalue; dx and ax are the zero-eigenfunctions 
(coUisional invariants) with lox — 0, normalized as {ax\dx) — 1. Here left and 
right eigenvectors are different because the matrix Q is asymmetric in lattice 
gas automata that violate the condition of detailed balance. The vectors ax 
and jx are related as jx — Qciax, where ci — k ■ c and Q projects out the 
components of ciax parallel to the coUisional invariants. In case the collision 
rules are purely local the matrix A vanishes. The sound damping constant F is 
a linear combination of kinematic viscosity v and bulk viscosity C given by, 

F = (l-i)^.+ iC 
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Figure 1: Damping constant zo{k) of diffusion mode versus k in negative diffu- 
sion model [ Alexander 1992| , where a is the angle between k and ci, the x-axis. 



Density p = 2.0 and threshold Tc — 3.8. At T = 4.0 (dashed lines) the modes 
are stable; at T = 3.5 (solid lines) the modes are unstable. 

In the standard detailed balance models the eigenvalues ujx are simple func- 
tions of the density. In models violating the Stueckelberg conditions, such as 



the FCHC lattice gas [Coevorden 1993 or the biased triangular lattice gas 



[ Bussemaker 1993a , the eigenvalues have to be calculated numerically. Al- 
though no H-thcorem can be proven for these models, it turns out that the 
eigenvalues still satisfy the inequalities, < < 2, although transport coeffi- 
cients may be negative (see section 3.2). 

Next we turn to the numerical solution of the eigenvalue equation (^). For 
a given set of collision rules we evaluate the matrix elements Q.ij{k) and the 
eigenvalues z\{k) numerically. Typical results are illustrated in Fig. 1 for a 
square lattice gas, described in section 3.3, with a single slow (diffusion) mode. 
The figure shows the damping constant Z]j{k) of the diffusion mode as a function 
of k. The diffusion mode is only stable (dashed lines) if the control parameter T 
(defined in section 3.3) exceeds a certain threshold value Tc- For T < Tc (solid 
lines in Fig. 1) the modes are unstable in two independent fc-directions. 

2.3 Cahn-Hilliard theory 

What is the important information that can be obtained from the spectrum 
zx{k) regarding stability of fluctuations and behavior of structure function and 
density-density correlation function? The criterion for dynamic stability of long 
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wavelength fluctuations is Kez\{k) < 0, which reduces to the criterion Fa > 
according to in cases where vx is reaL Furthermore, if there is a control 
parameter, say T, then the root Tc of the equation, Ilez\{k, T) = 0, determines 
the threshold value Tc, below/above which the A-th mode is unstable/stable. 

Numerical evaluation of the unstable eigenvalue z\{k) provides the following 
additional information: (i) the fc-direction of maximum instability (in Fig. 1: 
a = 45°); (ii) the cut-off wavelength Ac = 27r/fcc above which all modes are 
unstable (in Fig. 1: kc ~ 0.77 at a = 45°); (iii) the wavelength of maximum 
growth A„i = 27r/fc„j (in Fig. 1: fc,„ ~ 0.55 at a = 45°); (iv) the onset time 
of instability Tq — [Re z\{km)]~^ (in Fig. 1: Tq ~ 350 at a = 45°). The corre- 
sponding eigenfunctions are the order parameters, i.e. in shear modes it is the 
transverse momentum density and in sound modes it is a linear combination of 
longitudinal momentum and mass density. 

The existence of unstable long wavelength modes indicates that the system 
starts to phase separate and develop spatial structure on a typical initial length 
scale Am- The fc-directions of maximum growth determine the patterns in the 
structure function. Once eigenvalues and eigenfunctions have been computed, 
and possible unstable modes with Iiez\{k) > have been identified, one can 
use the Cahn-Hilliard theory to evaluate the structure function S{k, t) or its 
inverse Fourier transform, the density-density correlation function G{f,t), 

5(fc,i) = T/-i(|p(fc,t)|^) =^exp[-*fc .7^6(^0 (7) 

r 

where V is the number of nodes in the system, and p{k, t) the Fourier component 
of the microscopic density fluctuation 5p{r^ t) = p{r^ t) — p around the uniform 
density p = J^i fi- Its long wavelength components develop essentially like the 
hydrodynamic modes p{k, t) ~ e-icp[zx{k)t]. Hence the structure function at the 
onset of instability is determined by the unstable long wavelength modes, i.e. 

Q(t +\ { S{k,Q)ex^[2Rezx{k)t] [k < k^) 

^^''^'^ - \ Sik,0) (fc>fcc) ^ ^ 

The arguments presented here are essentially the Cahn-Hilliard theory of spin- 



odal decomposition | Langcr 1992 1 . We assume that the structure function S{k,t) 



= S{k,0) remains in equilibrium for k > kc- 

3 Phase separation 

3.1 Thermodynamic instabilities 

We consider a lattice gas model for a gas-liquid system with an attraction of 
(long) range R, which is the control parameter. It is essentially a 7-bit FHP 
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model to which a square well attraction of range R has been added. If the Fermi 
exclusion principle permits, there is an instantaneous exchange of momentum 
between two particles a distance R apart, such that their relative velocity is 
always reduced. The model therefore violates detailed balance. The momentum 
flux contains not only a kinetic part, but also a collisional transfer part. 

The perturbation theory of section 2 has been applied to this gas-liquid 



model iGerits 1993| . The pressure p, the speed of sound and the two vis- 
cosities and C have been calculated analytically, and have been succesfuUy 
compared with computer simulations. For the present discussion only the equa- 
tion of state, 

p = if-2Rf\l-f)\ (9) 

is of interest, as it shows the thermodynamic instability. Here 6f — p~ fo is the 
average number of moving particles per node, and fo that of rest particles. For 
R > Rc = 8 the equation of state exhibits a van der Waals loop. In the density 
interval (spinodal region) where dp/ dp — v1 < 0, the spatially uniform state 
is thermodynamically unstable. In computer simulations one observes droplet 
formation. The measured pressure p at short times {t < 500) is in agreement 
with (||) for all densities. For t > 500 one observes the coexistence pressure Pc- 

Inside the spinodal region (1.7 < p < 2.6) the pressure Pc is in excellent 
agreement with the pressure obtained from Maxwell's equal area construction. 
In the metastable regions, where nucleation is the basic mechanism for droplet 
formation, the coexistence pressure Pc could only be observed after artificially 
creating some high or low density droplets, acting as nuclei of condensation. 

The spinodal decomposition can also be analyzed in terms of the dynamic 
instabilities of section 2. Here the speed of sound is purely imaginary, Rc Za-{k) = 
±|a|fc. Hence the sound modes stop to be propagating, and the amplitude of 
one of them starts to grow as ~ exp[|a|fci]. Also in a continuous gas-liquid 
system spin odal decomposition is closely linked to the instability of sound modes 



[Koch 1982 



3.2 Unstable propagating modes 

The biased triangular ledtice gas, introduced by Bussemaker and Ernst [1993a], 
shows a dynamic instability with unstable sound modes that remain propagat- 
ing. The model is a stochastic 7-bit lattice gas with strictly local collision rules, 
that preserve the lattice symmetries, but violate detailed balance. Typically, 
the forward and backward collision probabilities are different whenever a rest 
particle is involved. 

The most characteristic feature of this model is shown in Fig. 2. It has been 
observed in computer simulations, and confirmed by numerical solution of the 
nonlinear lattice Boltzmann equation ([^) . There is a density interval, where the 
the spatially averaged density of rest particles fo{p) has a negative slope. In this 
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Figure 2: Density of rest particles fo{p) versus p in the biased triangular model. 



Mean field theory (solid line) and simulations (triangles) [Bussemaker 1993a . 



model the pressure is purely kinetic, 



P 



(10) 



and the speed of sound Vs, given by = ^[1 ~ dfo{p)/dp], can exceed the value 
1/V2 ~ 0.71. In the half-filled lattice Vs ~ 0.76, as can be read off from Fig. 2. 
Consequently the bulk viscosity in the biased triangular lattice gas is negative 
on account of (||). 

A negative bulk viscosity by itself is not a sufficient condition for instability, 
as can be concluded from the work of Bussemaker and Ernst [1992], where 
several stable lattice gases with the characteristics of Fig. 2 have been simulated. 
The sound modes only become unstable if the sound damping constant F = 
5(^ + becomes negative, i.e. if in addition v is sufficiently small. These 
conditions are realized in the biased triangular lattice gas. 

More quantitative information is obtained by solving the eigenvalue problem. 
One finds a density interval where sound modes are unstable. For instance, at 
p = 3.5 one has Ileza(k) > for k < kc with kc ~ 0.6; it reaches a maximum 
value = Reza{km) — 0.005 at km ~ 0.4 [Bussemaker 1993a]. 

The onset of phase separation, the coarsening phenomena, the late stages of 
growth, the scaling properties of the structure function S{k,t) for the different 
order parameters, and the exponent a of the dynamic correlati on length Am ~ ^" 
in this model, have been extensively discussed in the literature | Bussemaker 1993q , 
Bussemaker 1993bl. 
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3.3 Fluctuating bias fields 



The best known model of this kind is the Rothman-KeUer model [ Rothman 1988 | 
of a binary fluid consisting of red and yellow particles {a = r,y), which are me- 
chanically identical. There is an attraction between like particles to enhance 
phase separation. The basic model is again a triangular lattice gas with seven 
different velocity channels and at most one particle per channel, i.e. the occu- 
pation number Si{r) = ot 1. To indicate the color {a — r, y) one adds a color 
label to the occupation number s^i (r) with Sri {r) -\- Syi (r) = Si (r) . 

The transition probability A{s s') from a pre-collision configuration 
s(r) = {scri(r); a = r,y;i = 0,1,2, ... ,b} at node r to a post-collision configu- 
ration s'{r) contains the usual delta functions accounting for the conservation 
per node of red and yellow particles and of total momentum. In addition there 
is a bias factor exp[/?J(s') • G(s„.„.)]. The normalization is fixed by the con- 
dition J2s' ~^ — 1- The bias factor represents the effect of a local bias 
field G(s„.„.), which is here a color gradient, that depends on the configurations 
{Sn.n.} of all nearest neighbor nodes. It further contains the color flux J(s') in 
the post-collision configuration s'(r). 

The quantity AW = J(s') • G(s„.„.) represents the amount of work done 
on the system by sending particles up a color gradient. The control parameter 
T = is a temperature-like variable. The relative probability for the transition 
s ^ s' is given by the bias factor exp[/3Aiy]. Depending on the value of /3AW 
transitions are less or more likely. At very high temperature (T — > oo) there 
is no bias; at very low temperature (T 0) there is a strong bias. There is a 
threshold value Tc below which the color diffusion mode is unstable. The order 
parameter is the difference in concentration of red and yellow particles. 

Similar fluctuating bias fields are present in the negative viscosity model of 
Rothman [1989]. The fluctuating bias field is the gradient of the local flow fleld 
and the flux is the momentum flux. This model leads to spontaneous ordering 
in the velocity field with regions of high vorticity. The order parameter is the 
transverse momentum density or vorticity. 

A further simplification of the Rothman-KeUer model is the negative dif- 
fusion model Alexander 1992| . It is defined on a square lattice. Total mass 



per node is conserved, total momentum is not. The model has only a single 
slow mode. The bias factor exp[/3J(s') • G(s„.„.)] contains the gradient of the 
microscopic density p(s(r^) and the particle current, 

G(s„.„.) = y^CmP{s{r + c„)) 

rn 

/(.') = (11) 

i 

The collection of models described above violate the condition of detailed bal- 
ance. 
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Figure 3: Negative diffusion model at density p — 2 and temperature T = 3.5. 
No external driving field, (a) Unstable (positive) parts of zi:,{k) around origin 
of Brillouin zone, denoted by shades of gray. Stable (negative) values are left 
blank. The spectrum has the symmetry of the square lattice, (b) Density- 
density correlation function G{r, t) ait — 100, obtained from the eigenvalues in 
(a) using Cahn-Hilliard theory. Maximal positive values are denoted by black, 
maximal negative values by white, and intermediate values by shades of gray, (c) 
Simulated values of G{r, t) aX t = 100, obtained from one run for a 1024 x 1024 
system at T = 3.0. 



We need to caution the reader here, because it is not at all clear whether 
the forces driving the kinetics of phase separation in diffusive systems or the 
kinetics of ordering of the flow field in regions of high vorticity, can be sensibly 
modeled by sending fluxes upstream against the prevailing gradients (negative 
transport coefficients). According to irreversible thermodynamics fluxes are 
downstream, transport coefficients are positive, as is the irreversible entropy 
production. These questions need further clarification. 

Here we will only illustrate some of the dramatic effects, such as pattern for- 
mation, that results from negative transport coefficients induced by fluctuating 
bias fields. 

As an illustration we discuss the negative diffusion model and apply the 
linear stability analysis of section 2. For the special case of the half -filled lattice 
(yO = 2) the eigenvalue equation (^ can be solved analytically to yield the 
dispersion relation, 

Z£)(fc) — log[i(cos /cj; -f coaky) + 4ci;(sin^ + sin^ ky)] (12) 

for the diffusion mode, where w is a typical matrix element of the nonlocal 
collision operator, defined by Alexander et al. [1992]. We first apply the stability 
criterion of section 2. 

The eigenvalue Z]j{k) is plotted in Fig. 1, for values of the control parameter 
above (T — 4.0) and below (T — 3.5) threshold. In Fig. 3a the unstable regions 
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of ZD{k) in the (kx, A:y)-plane are indicated in shades of gray. The maxima are 
located in the 45°-direction, a typical distance km away from the origin. As T 
is further decreased below threshold, the typical wavelength Am = '^t^ jkm of the 
most unstable excitation decreases. The threshold value follows directly from 
the small-A; expansion of (p^), yielding Z]j{k) = — (i — 4(jj)fc^ = —Dk^. For 
w > Wc = 1/16 (which corresponds to T < ~ 3.8 according to Bussemaker 
and Ernst [1993c]) the diffusion coefficient takes a negative value, and the system 
becomes unstable for long wavelength fluctuations. The mean field value found 
for the threshold (Tc ~ 3.8) is about 20% higher than the threshold measured 
in the computer simulations of Alexander et al. [19921. 

According to the Cahn-Hilliard theory we can use (0) and (H) to calculate the 
initial structure and patterns in the density-density correlation function G'(r, t). 
This yields the checkerboard pattern of Fig. 3b with a typical lattice distance 
Am. Wc have also carried out computer simulations of the structure function 
and its inverse G{r, t) for a 1024 x 1024 system prepared in a random initial state 
at t = 0. Quantitative comparison with the Cahn-Hilliard prediction of G(f,t) 
is complicated by the 20% discrepancy in the threshold temperature, mentioned 
in the preceding paragraph. Fig. 3c shows G{f,t) at t = 100, obtained from a 
simulation at T = 3.0, a value which is just below the simulation value of the 
threshold temperature Tc =3.05-3.10 reported by Alexander et al. [1992]. One 
sees that the structure of G(r, t) in the actual simulations is qualitatively the 
same as in the Cahn-Hillard theory, although the typical length scale is about 
30% smaller. 



4 Striped phases 

4.1 Driven diffusive systems 

If the spatial symmetry of the underlying lattice is broken by an external driving 
field or by a spatially uniform flow, the /c-regions of the most unstable excitations 
become strongly anisotropic, and different types of striped phases may appear. 
The occurrence of striped phases in driven diffusive systems is well known 



[Zia 1993 1 . They have also been observed in computer simulations of the nega- 



tive diffusion model, where an external field F is added, that drives a particle 



current [Alexander 1992 1. The driving force in lattice gases can be implemented 
by replacing the bias factor in section 3.3 by exp[ (/3G(s„.„.) -I- F) ■ J(s') ]. 

In our analysis we restrict ourselves to the linear response regime, where the 
particle current satisfies the constitutive relation, 

J^fiF + DVp (13) 

with fj, the mobility and D the diffusion coefficient. Both transport coefficients 
are related by the Einstein relation D = Tfi. 
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Figure 4: Negative diffusion model at density p = 2, in the presence of an 
external driving field F — (0.5; 0.0). (a) Unstable (positive) parts of ZD(k) 
around origin of Brillouin zone, denoted by shades of gray, at T = 3.6. Stable 
(negative) values are left blank. The symmetry of the square lattice is broken 
by the driving field, (b) Density-density correlation function G(r, t) ait = 100, 
obtained from the eigenvalues in (a) using Cahn-Hilliard theory, (c) AtT — 3.5 
the eigenvalue zoik) has become unstable in all fc-directions. 



For sufficiently high temperature T there exists a stable spatially uniform 
state. Its distribution function f°{F) can be obtained by numerically solving 
(0). The velocity distribution is anisotropic due to the presence of the field. As T 
is decreased the density fluctuations of long wavelength become again unstable. 
The stability analysis of section 2 is performed by linearizing (0) around f° (F) , 
and determining when and where the eigenvalue zoik) becomes negative. As 
the spatial symmetry is broken the threshold of stability Tc{F, k) depends not 
only on the field F but also on the direction of k. 

A systematic analysis of the fc-regions of instability as a function of T will be 
presented elsewhere [ Bussemaker 1993c | . Here we only illustrate the occurrence 
of striped patterns for a typical case with F in the x-direction. At T = 3.6 
the eigenvalue zoik) only shows unstable modes if k is perpendicular to F, 
with peaks at \km\ — 0.2. Directions parallel to F are stable. See Fig. 4a. 
By combining the numerical results for Z]j{k) with the Cahn-Hilliard theory 
in (^ and (||) one obtains the density-density correlation function of Fig. 4b, 
indicating a striped pattern. 

Alexander et al. [1992] mention that "there also appears to be a second 
transition at a lower temperature to a phase with structures resembling those 
found for = 0" . Our analysis suggests that this transition is very gradual: as 
the temperature decreases further below Tc, the region of unstable fc-directions 
spreads out, until at some temperature all directions have become unstable 
(Fig. 4c). Of course for very low temperatures the influence of F becomes 
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negligible, since the term (3G then dominates the bias factor. 



4.2 Unstable uniform flows 

In computer simulations on FCHC lattice gases Henon [1992] has observed that 
an initially uniform flow is unstable. The flow orders itself into parallel stripes 
II uq with alternating flow velocities parallel and anti-parallel to Uo- The lattice 
gas violates the Stueckelberg condition of semi-detailed balance. The above ob- 
servations suggest that uniform flows tend to destabilize finite speed equilibria. 

To obtain some indications about possible instabilities in uniform flows in 
lattic e gases we si mplify the nonlinear collision operator in (]^) to a BGK-coUision 
term IChen 1992|| , i.e. 

= (14) 

T 

where all non- vanishing eigenvalues ujx — 1/t are related to a single relaxation 
time r. The 'local equilibrium' distribution function in a rf-dimensional lattice 
with lattice distance \ci \ = Vg has been consistently chosen as, 

f! - 




fo - /o-M^ (15) 



where p = fo + bf . This choice guarantees that the momentum flux in local 



equilibrium is Galilei-invariant up to and including C'(w^)-terms, i.e. 



P^a,i3 + pUaUjj where a and (3 denote Cartesian components and p = (w^ / d) bf is 
the hydrostatic pressure, as in (p^. This choice still leaves fo as a free parameter 
to model the speed of sound through = {v'^/d)[l — dfo{p) / dp]. Alternatively 
one may add extra collision terms on the right hand side of ( p^ to model the 
mechanism of collisional transfer. 

In either case one can apply the method of section 2 to investigate the 
stability of the uniform stationary state ff{uo) with a flnite flow velocity Uq- It 
turns out that the shear modes remain stable, but the sound modes (cr — ±) 
may become unstable. The sound damping constant is calculated as 

r.(fc) = i^ir - i) + auoe) { 4^ - (v, + aUoeA (16) 



2vs" ' {d + 2 

where Uoi — k ■ Uq. The shaded area in Fig. 5 gives the domain in the (ws, Uot)- 
plane where both sound modes are stable. The lines marked {a = +) and 
(cr = — ) denote the zero's of (p^. The A;-direction parallel to Uq is the most 
unstable [|. 

^ After completion of this research we received the preprint on ' StabiUty analysis of lattice 
Boltzmann methods' by J.D. Sterling and S. Chen, in which stability diagrams similar to 
Fig. 5 have been obtained by numerical solution of BGK-models. 



12 




Figure 5: Stability diagram for uniform flow in {vg, Uo£)-plane. In the shaded 
region both sound modes are stable. 



If Uoi = we recover the BGK-approximation to the results (|6|), because 
all non-vanishing eigenvalues satisfy uju = ijJc = 1/t. If Vs > Vo[3/{d + 2)]^/^ 
the sound modes are unstable in all directions, and the BGK-lattice Boltzmann 
equation constitutes a mathematical model for the unstable propagating sound 
modes, discussed in section 3.2. 

Inspection of the diagram in Fig. 5 shows that it is easy to construct BGK- 
models where sound waves with k \\ Uo are unstable, and those with k Jl Uo are 
stable. In such models the dynamic instability leads again to phase separation 
and striped patterns, oriented perpendicular to the direction of the flow. The 
order parameter is the sound mode, which is a linear combination of momentum 
density parallel to Uo and mass density. 

The diagram of Fig. 5 indicates how a stable basic equilibrium {uo — 0) can 
become unstable when given a non- vanishing total momentum P — Nuo, and 
how striped patterns can be formed with the sound mode as order parameter. 
However, we have not succeeded in constructing a BGK-model, that shows the 
Henon instability. Adding a flow velocity Uq to the negative viscosity model of 
section 3.3 may lead to an unstable shear mode at sufficiently high flow velocity, 
but generates stripes perpendicular to the flow. 

In summary^ we conclude that the stability analysis based on mean field theory 
gives qualitative and quantitative information about the initial patterns and 
accompanying length and time scales in phase separation problems. The method 
presented here is applicable to both the Boltzmann equation for lattice gas 
automata, and the mathematical model Boltzmann equations, referred to as 
BGK-models. 
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